本发明属于生物信息技术领域,尤其涉及一种高通量转录组测序数据的分析方法。所述分析方法更加有效准确的分析样本当前的基因转录情况,其在分析过程中通过Burrows-Wheeler变换方法定位后,再通过Hash表进行了二次定位。本发明通过WGA加Hash alignment的方法给出了一种更加可靠有效的高通量转录组测序数据的分析方法。
1.一种高通量转录组测序数据的分析方法,其特征在于,其在分析过程中通过Burrows-Wheeler变换方法定位后,再通过Hash表进行了二次定位,所述分析方法的具体步骤如下:(1)WGAxa0xa0Alignment方法:获得RNA-Seq的原始数据后,将转录组测序数据以可变长度和步长切成短片段,并应用Burrows-Wheeler变换方法定位到参考基因组序列上,匹配允许1个碱基的错配,并设定最大允许查找定位的数目;所述可变长度与步长可以根据实际Reads的长度来确定;(2)将匹配上的测序数据分为节点和数据链,基于参考基因组mRNA注释信息以及匹配度打分来滤除不可靠节点和数据链;(3)切割所有匹配的参考基因组序列构成新的有效的重构参考序列,应用hash表进行二次定位,以此提高匹配准确度及速度;(4)对Reads匹配结果重新定位到原始参考基因组序列上,重复步骤(2),确定Nodes、link以及Chain,基于基因组注释mRNA信息以及数据链各exon间连接度找到新的mRNA可变剪切、基因融合和新基因,并给出基因表达水平、变异情况的统计信息;(5)基因Chain和Nodes进行分层展示,并根据不同Node、link类型给出不同颜色分类标识;(6)生成各类基因数据报表,并建立索引,实现快速定位。
2.根据权利要求1所述的分析方法,其特征在于,所述步骤(2)中,具体步骤如下:(a)Node确定:将基因组上有Reads覆盖的连续区域定为一个Node,根据设定的基覆盖Coverage;(b)link确定:I.载入所有形成link的reads,也就是指单条Read定位到基因组的不连续的两段以上区域;II.进行峰值处理,依据聚类分析结果给出link的起始位置,并根据如下规则进行修正:当Reads覆盖度大于10,以超过50%形成link的Reads的起始位置进行修正;当Reads覆盖度小于10,以超过70%形成link的Reads的起始位置进行修正;III.依据进行基因注释信息和GT-AT法则进行link起始位置修正,去除边缘可变剪接Nodes的影响,对边缘1-3bp碱基进行exon与exon剪接接头检查,去除重复碱基定位的误差;IV.小Gapxa0link的检测:设定最小exon间碱基数目为6,低于6bp的Gap进行删除变异处理,而大于6bp的Gap定义为Gapxa0link;(c)Chain确定:连续的Node及link定义为一个Chain;I.根据Node及连接Node的link,给出所有可能的不重复的基因Chain;II.根据基因注释信息,找出基因区域内的所有可能的基因Chain;找出基因区域内的所有的新基因Chain;III.滤除不可靠基因Chain:相似基因Chain、重叠基因Chain以及融合基因chain按各自规则依次进行滤除;相似基因Chain滤除规则:序列相似度大于80%的基因进行相似基因Chain检查,将相似基因间的相同基因Chain进行合并,保留不同的基因Chain;重叠基因Chain滤除规则:如果重叠部分和某一基因的mRNA完全匹配则直接划分到这一基因,如果没有完全匹配的,则比较所在chain和各个基因的mRNA匹配的长度,基因Chain定位到匹配长度最长的基因;融合基因Chain滤除规则:对于相邻近基因只保留相同模板链上相同方向的基因融合Chain;重叠基因间的融合基因Chain滤除,根据覆盖度拆分定位到对应的基因上;滤除相同染色体上覆盖范围大于10000的融合基因Chain;支持多基因间融合查找。
3.根据权利要求1所述的分析方法,其特征在于,所述步骤(3)中,具体步骤如下:(a)xa0xa0xa0xa0以所有基因组注释基因的mRNA,有Reads覆盖的生成基因chain的序列为基本参考序列;(b)xa0xa0xa0设定提取序列初始长度为100,再动态计算平均Reads长度值后替换;(c)xa0xa0xa0xa0从基本参考序列中分别取exon与exon的接头序列、基因Chain的Node与Node间的连接序列,长度等于提取序列长度;去重复追加到重构参考基因组序列,并记录在原始基因组中位置以备还原定位;(d)xa0xa0xa0从基本参考序列中分别取exon中间序列、基因Chain的Node中间序列,去重复存为重构参考基因组序列,并记录在原始基因组中的位置以备还原定位;(e)xa0xa0xa0xa0拼接序列时减少重复序段的干扰,根据exon长度以及Node长度调整提取序列的大小;(f)xa0xa0xa0xa0在重构参考基因组序列中进行去重处理,以降低参考序列大小,提高运行速度;对测序样本用Hashxa0Alignment方法对重构参考基因组序列进行匹配,可设定最小匹配长度以及最低相似度,默认设置为12及85%。技术领域
本发明属于生物信息技术领域,尤其涉及一种高通量转录组测序数据的分析方法。
背景技术
随着新一代高通量测序技术的快速发展,建立在高通量测序基础上的转录组测序技术已成为目前从全基因组水平研究基因表达和转录组分析的重要手段.转录水平的调控是生物体最主要的调控方式.在深度测序技术出现之前,高通量测定不同基因转录水平的主要手段是基因芯片,它可以对不同组织或不同发育阶段的基因表达差异和模式进行分析,而RNA-Seq技术最基本的应用也是检测基因的表达水平,它对同一样品深度测序可以捕获低表达的基因,而对大量样品同时测序可以获得样品之间的表达差异。与基因芯片数据比较,RNA测序得到的是数字化的表达信号,无需设计探针,能在全基因组范围内以单碱基分辨率检测和量化转录片段,具有灵敏度高、分辨率高和应用范围广等优势。除此之外,xa0研究人员还可以获得转录本表达丰度、转录起始位点和可变剪切等重要信息。所以,建立在高通量测序基础上的转录组研究已经逐步取代基因芯片技术成为目前从全基因组水平研究基因表达的主流方法。在对基因组进行测序时,基因组被切割成很多小片段,通过复制、碱基辨识等步骤,我们可以获得这些短序列的碱基序列(测序序列)。然而在切割基因组后,我们无法知道各个测序序列的相对位置。如果没有参考基因组,就只能通过装配技术来得到所测的基因组。如果有一个已被测得的基因组作为参照,这就是一个相对容易的重测序问题。现在我们在生物学研究、个体化医疗中面临的测序问题,绝大部分是或可以近似转化为重测序问题。在重测序问题中,我们要寻找每一个测序序列在参考基因组上的位置或坐标,我们称之为测序序列定位。例如,通过将转录组的测序序列定位至参考基因组序列来测量不同种类的RNA的含量;或者通过将转录组的测序序列定位至参考基因组序列来探测mRNA可变剪切的模式等。通常的高通量测序数据的定位(Mapping)方法,主要有基于哈氏表(hash)的方法和基于Burrows-Wheeler变换的方法。基于哈氏表(hash)的方法只能处理定长的序列数据,适用于参考序列较短的;如果序列太长(如人的全基因组序列),则需要很大的计算机内存,处理速度也很慢。优点是较好地支持容错率。基于Burrows-Wheeler变换的方法能处理变长的序列数据,但不能支持较大的容错率,否则将丧失速度优势。由于测序仪器自身的测序误差,实际的高通量测序数据往往有很多的噪声,这就给正确的定位带来很大的困扰。对于真核基因组情况将更加复杂,由于真核基因是断裂基因,由内含子和外显子构成,转录机制更加复杂。mRNA是按照特定的转录机制由外显子序列拼接而成,也就是说一条转录组测序数据(Read)定位在基因组参考序上可能是不连续的片段,这就为正确定位测序数据、基因转录水平、可变剪接以及基因融合带来更大的难度。
发明内容
本发明目的在于针对现有技术的不足,提供一种高通量转录组测序数据的分析方法,更加有效准确的分析样本当前的基因转录情况,其在分析过程中通过Burrows-Wheeler变换方法定位后,再通过Hash表进行了二次定位。在本发明的一个实施方案中,所述分析方法的具体步骤如下,实施流程图见附图1:(1)WGAxa0xa0Alignment方法:获得RNA-Seq的原始数据后,将转录组测序数据以可变长度和步长切成短片段,将转录组测序数据以可变长度和步长切成短片段,并应用Burrows-Wheeler变换方法定位到参考基因组序列上,匹配允许1个碱基的错配,并设定最大允许查找定位的数目,设定有效值1-50;可变长度默认值为30,步长默认值为5,所述可变长度与步长可以根据实际Reads的长度来确定,Reads长度越长,可变长度及步长可以设定得越大,以提高运行效率,步长设定为可变长度的1/6为佳;(2)对匹配上的测序数据分为节点和数据链,基于参考基因组mRNA注释信息以及匹配度打分进行滤除不可靠节点和数据链,具体步骤如下:(a)Node确定:将基因组上有Reads覆盖的连续区域定为一个Node,根据设定的覆盖深度保留符合条件的Node,可分别已知基因区最小Read覆盖深度和非已知基因区最小Read覆盖深度,默认值分别为1和4。根据不同测序数据质量可以进行调整,测序深度越大,设定得最小Read覆盖深度以滤除噪声。(b)link确定:I.载入所有形成link的reads,也就是指单条Read分别定位到基因组的不连续两段以上区域;II.进行峰值处理,依据聚类分析结果给出link的起始位置(参考文献1),并根据如下规则进行修正:当Reads覆盖度大于10,以超过50%形成link的Reads的起始位置进行修正;当Reads覆盖度小于10,以超过70%形成link的Reads的起始位置进行修正(实施案例见附图12);III.依据进行基因注释信息和GT-AG法则进行link起始位置修正,去除边缘可变剪接Nodes的影响,对边缘1-3bp碱基进行exon与exon剪接接头检查,去除重复碱基定位的误差;VI.小Gapxa0link的检测:设定最小exon间碱基数目为6,低于6bp的Gap最为删除变异处理,而大于6bp的Gap定义为Gapxa0link。(c)Chain确定:连续的Node及link定义为一个Chain。I.根据Node及连接Node的link,给出所有可能的不重复的基因Chain。II.根据基因注释信息,找出基因区域内的所有可能的基因Chain。找出基因间区域的所有的新基因Chain。III.滤除不可靠基因Chain:相似基因Chain、重叠基因Chain以及融合基因chain(该基因Chain覆盖区域跨越两个以上基因的区域)按各自规则依次进行滤除。相似基因Chain滤除规则:序列相似度大于80%的基因进行相似基因Chain检查,将相似基因间的相同基因Chain进行合并,保留不同的基因Chain。重叠基因Chain滤除规则:如果重叠部分和某一基因的mRNA完全匹配xa0则直接划分到这一基因xa0如果没有完全匹配的xa0则比较所在chain和各个基因的mRNA匹配的长度,基因Chain定位到匹配长度最长的基因。融合基因Chain滤除规则:对于相邻近基因只保留相同模板链上相同方向的基因融合Chain;重叠基因间的融合基因Chain滤除,根据覆盖度拆分定位到对应的基因上;滤除相同染色体上覆盖范围大于10000的融合基因Chain;支持多基因间融合查找。(3)切割所有匹配的参考基因组序列构成新的有效的重构参考序列,应用hash表进行二次定位,以此提高匹配准确度及速度;拼接参考序的步骤如下,实施见附图2:(a)xa0xa0xa0xa0以所有基因组注释基因的mRNA,有Reads覆盖的生成基因chain的序列为基本参考序列;(b)xa0xa0xa0设定提取序列初始长度为100,再动态计算平均Reads长度值后替换;(c)xa0xa0xa0xa0从基本参考序列中分别取exon与exon的接头序列、基因Chain的Node与Node间的连接序列,长度等于提取序列长度;去重复追加到重构参考基因组序列,并记录在原始基因组中位置以备还原定位;(d)xa0xa0xa0从基本参考序列中分别取exon中间序列、基因Chain的Node中间序列,去重复存为重构参考基因组序列,并记录在原始基因组中的位置以备还原定位;(e)xa0xa0xa0xa0拼接序列时尽量减少重复序段的干扰,根据exon长度以及Node长度调整提取序列的大小;(f)xa0xa0xa0xa0在重构参考基因组序列中进行去重处理,以降低参考序列大小,提高运行速度;(g)xa0xa0xa0对测序样本用Hashxa0Alignment方法对重构参考基因组序列进行匹配,可设定最小匹配长度以及最低相似度,默认设置为12及85%;(4)对Reads匹配结果重新定位到原始参考基因组序列上,重复步骤(2),确定Nodes、link以及Chain,基于基因组注释mRNA信息以及数据链各exon间连接度找到新的mRNA可变剪切,基因融合和新基因,并给出基因表达水平、变异情况的统计信息;(5)xa0基因Chain和Nodes进行分层展示,并根据不同Node、link类型给出不同颜色分类标识,见附图3、附图5-12,具体规则如下:分层展示原则:与已知mRNA吻合的基因显示在第一层;发生可变剪辑变异的Node、独立的新link放在第二层;重叠基因区域的按照前两个原则顺序分层显示。颜色分类标识定义如下:已知基因Nodes:与已知基因exon完全吻合。紫色横线标识,见附图3。已知基因link:与已知基因的mRNA剪接方式相吻合的link。用紫色连接线标识,连线上方数字为形成该link的Reads数目,见附图3。新Node(Novelxa0Node):该Node与任何已知注释基因的exon区域没有重叠区域,蓝色标识。新link(Novelxa0link):连接link的至少一个Node是新Node;连接link的至少一个Node的连接位置与已知基因exon位置有差异。用蓝色连线标识,连线上方数字为形成该link的Reads数目,见附图5-6。融合link(Fusionxa0link):两个不同基因见发生可变剪接的link,不直接标识,只在数据报表中给出link的起始位置,点击报表可以在基因Chain的显示图中给出定位区域。见附图7。可变剪接(Alternativexa0Splicing):与已知mRNA的exon不完全重合的,至少有一端与exon的位置有差异,用桔黄色连线标识,见附图5,附图8。可变转录起始(Alternativexa0Transcriptxa0Start):与mRNA第一个exon不完全重合的Node,用灰色横线标识,见附图9可变转录终止(Alternativexa0Transcriptxa0Stop):与mRNA最后一个exon不完全重合的Node,用灰色横线标识,见附图6,附图10插入Node(Insertion):在基因内与已知exon没有任何重叠的Node,并且没有完全跨越相邻的两个exon,用红色横线标识,见附图8。通读(Readxa0Through):跨越两相邻基因间区域的Node,用灰色横线标识,见附图11。内含子包容(Intronxa0Retention):扩越基因mRNA两个已知相邻exon的Node,用红色横线标识,见附图12。非mRNA的Node:除了mRNA之外的RNA的Node,包括miscRNA,xa0ncRNA等,用绿色横线标识,见附图6。(6)生成各类基因数据报表,并建立索引,实现快速定位,见附图3和附图4展示。数据报表包括每个Nodes、link的染色体定位信息、基因信息、Reads覆盖深度、形成link的Reads数目、Pair-Endxa0Reads通过形成link的数目、已近Node的类型、Isoform和蛋白质的NCBI识别号,可直接链接到NCBI网站查看mRNA和蛋白质的详细信息。数据报表(附图4)与基因Chain展示(附图3)实现关联定位,可直接定位到需要查看的基因、Nodes、link以及Chain,也可根据基因名字、染色体位置和基因序列直接进行搜索。
附图说明
图1:本发明高通量转录组测序数据的分析方法的流程图图2:本发明重构参考序方法示意图图3:本发明转录组测序数据分析结果显示图,C部分为基因Chain显示,其中A为已知Node(紫色),B为已知link(紫色),该图显示的为获得的一条已知基因Chain(紫色),D部分为基因注释信息显示,其中G为基因区域,H为exon,该基因有6个mRNA,E部分为参考基因组的整体显示,双击可定位到相应基因区域,F部分为Reads覆盖深度显示。图4:本发明转录组测序数据分析结果报告图图5:xa0本发明发现新link(Novelxa0link)及可变剪接(Alternativexa0Splicing)示例,图中三个标注为A的Node为可变剪接(桔黄色),标注为B的link为新link(蓝色)。标注为C的为已知Node(紫色),标注为D的为已知link(紫色)图6:本发明发现新link及MiscRNA示例,图中标注为A的三个Node是MiscRNA的Node(绿色),标注为B的link是新link(蓝色),标注为C的link是已知link(紫色),标注为D的Node是可变转录终止Node(灰色)。图7:本发明发现融合link(Fusionxa0link)示例,图中A区域中大斜线的区域为融合link的定位,B区域第一条记录给出了融合link发生的起始位置(Start,xa0End)、染色体(Chr)、基因(Gene)、形成该link的Reads数目(linkxa0Number)等信息。图8:本发明发现插入Node(Insertion)示例,图中标注为A的Node为插入Node(红色),标注为B的Node为可变剪接(桔黄色),因为两个Node中间是不连续的。标注为C的两个link为新link(蓝色),标注为D的link为已知link(紫色),标注为E的Node为已知Node(紫色)。图9:本发明发现可变转录起始(Alternativexa0Transcriptxa0Start)示例,图中标注为A的Node为可变转录起始(灰色),标注为B的link为已知link(紫色),标注为C的Node为已知Node(紫色)。图10:本发明发现可变转录终止(Alternativexa0Transcriptxa0Stop)示例,图中标注为A的Node为可变转录终止(灰色),标注为B的link为已知link(紫色),标注为C的Node为已知Node(紫色)。图11:本发明发现通读(Readxa0Through)示例,图中标注为A的Node为通读(灰色),标注为B的Node是YOS1基因的可变剪接(桔黄色),标注为C的Node是PTP3基因的已知Node(紫色)。图12:本发明发现内含子包容(Intronxa0Retention)示例以及滤噪示例,图中标注为A的Node是内含子包容(红色),标注为B的Node为可变剪接(桔黄色),标注为C的Node为已知Node,标注为D的为已知link。依据本文步骤2中对link起止位置的进行了修正,给出了正确的link,并给出了同时存在的可变剪接及内含子包容。
具体实施方式
下面将结合附图以及进一步的详细说明来举例说明本发明。需要指出的是,以下说明仅仅是对本发明要求保护的技术方案的举例说明,并非对这些技术方案的任何限制。本发明的保护范围以所附权利要求书记载的内容为准。术语解释在本发明中所述“Read”是指高通量测序得到的一条序列。在本发明中所述“Node”是指reads定位到参考序后由连续的reads构成的结点。在本发明中所述“link”是指reads覆盖到两个以上node形成的连接。在本发明中所述“Chain”是指由nodes及经过nodes的link形成的链。在本发明中所述“Reads覆盖度”是指定位到参考序的连续区域上Read的数目。实施例1为使本发明的目的、技术方案和优点更为清楚,本发明用一套实际的人的高通量转录组测序数据为例进行数据分析:1.xa0xa0xa0xa0xa0xa0样本数据Illumina20由Illummina测序仪完成测序,转成fasta数据格式大小为3.29G,平均测序深度为30。2.xa0xa0xa0xa0xa0xa0以人的GRCh37.p5基因组序列作为参考序构建Index,序列文件大小为2.73G。3.xa0xa0xa0xa0xa0xa0依照表1设置参数,依照发明内容中的步骤1,以WGA比对测试样本Illumina20与人的参考基因组序列,匹配得到样本序列在参考序中的最佳匹配位置。4.xa0xa0xa0xa0xa0xa0依据发明内容中的步骤2的方法,确定所有Node,link以及Chain,去除不可靠link以及不可靠Chain。结果在表2中给出。5.xa0xa0xa0xa0xa0xa0依据发明内容中的步骤3描述方法,根据获得可靠Chain,拼接得到Hashxa0alignmentxa0方法的参考序,序列文件大小为85.5M。6.xa0xa0xa0xa0xa0xa0依据发明内容中的步骤4描述方法,对测试样本对新拼接得到的参考序进行Hashxa0alignment,比对参数见表1.7.xa0xa0xa0xa0xa0xa0依据发明内容中的步骤5描述方法,将比对结果重新定位到人的原始参考基因组序列上。8.xa0xa0xa0xa0xa0xa0再次依据发明内容中的步骤2的方法,确定所有Node,link以及Chain,去除不可靠link以及不可靠Chain。结果在表2中给出。9.xa0xa0xa0xa0xa0xa0基于基因组注释mRNA信息以及数据链各exon间连接度找到新的mRNA可变剪切,基因融合和新基因。10.xa0xa0依据发明内容中的步骤6描述方法,用不同颜色标识各类Node,link以及Chain,进行分层展示,快速查找及定位。11.xa0xa0统计基因表达情况,变异情况,以及新基因的情况。12.xa0xa0结果显示见表213.xa0xa0结果评价:尽管WGA方法可以快速实现测试样本序列的快速定位,仅用了66分钟,也能找到可靠的Node,linkxa0以及Chain,但根据与Hashxa0alignment结果比较可以明显看出,经过重新构建的基因组参考序,找到的Node总数较WGA方法减少了近1/3,但却可以获得相近的link数,说明并没有因为参考序列的大幅减少而丢失可靠的link。同时找到基因Chain不到WGA方法的一半,而获得的表达基因数目与实验情况基本吻合,同时融合基因定位比较合理,可见WGA方法获得的基因Chain受Nodes数目较大影响,有很大的孤立的基因Chain存在,而找到的47,627基因虽然没有给出详细划分,但人的基因总数目在3万左右,可以推测其中多数为孤立Node的或者单一link都成Chain的新基因,可靠性不高。因此可以清楚地确定尽管Hash方法以牺牲适度的时间换取了基因表达情况分析的可靠度,本发明通过WGA加Hashxa0alignment的方法给出了一种更加可靠有效的高通量转录组测序数据的分析方法。表1xa0参数设置ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="415"/>表2xa0WGA及Hash分析结果比较ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="411"/>参考文献1xa0肖宇,于剑,基于近邻传播算法的半监督聚类,软件学报,2008,19(11):2803-2813本发明内容仅仅举例说明了要求保护的一些具体实施方案,其中一个或更多个技术方案中所记载的技术特征可以与任意的一个或多个技术方案相组合,这些经组合而得到的技术方案也在本申请保护范围内,就像这些经组合而得到的技术方案已经在本发明公开内容中具体记载一样。